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A 2-D version of the asymmetric exclusion model for granular sheared flows is presented. 
The velocity profile exhibits two qualitatively different behaviors, dependent on control 
parameters. For low friction, the velocity profile follows an exponential decay while for large 
friction the profile is more accurately represented by a Gaussian law. The phase transition 
occurring between these two behavior is identified by the appearance of correlations in the 
cluster size distribution. Finally, a mean-field theory gives qualitative and quantitative 
good agreement with the numerical results. 
PACS numbers: 45.70.M, 05.10.G and 64.60.C. 

I. INTRODUCTION 

Among the numerous problems dealing with granular materials, one of the most challeng- 
ing is granular sheared flow Thus, although granular materials might exhibit solid, 
liquid or gas-like properties the granular flows cannot be described simply. Granular 
sheared flow arises in many different contexts such as pipe flow, pyroclastic flows , or even 
traffic jams Experiments performed in 2-D and 3-D geometries show surprising velocity 
profiles. First, the flow occurs only in a sheared region whose width is typically on the order 
of ten particle sizes (the shear zone). Also, velocity profiles appear to behave differently in 
2 and 3 spatial dimensions. The velocity decreases exponentially with the distance to the 
wall in 2-D with a small Gaussian correction [^], while in 3-D the profile is almost purely 
Gaussian 0. The goal of this paper is to present a two dimensional "toy model" where the 
velocity profile evolves correspondingly from an exponential like form to an almost Gaussian 
form. The model consists of vertically coupled layers. Each layer follows the well known 
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asymmetric exclusion (ASEP) model. In one dimension, the ASEP model has been widely 
studied and under certain conditions, exact solutions have been found using the infinite di- 
mension matrix method PJ. However, to our knowledge, very little is known concerning 2-D 
ASEP model. Our numerical simulations will in fact show a cross-over between an expo- 
nential velocity profile and a Gaussian velocity profiles when the control parameter crosses 
the value 1/2. This could be interpreted as a phase transition in infinite size systems, as 
the study of the clusters size distribution will indicate. Finally a mean field approach will 
be developed for the low values of the control parameter. 

II. THE MODEL 

The 1-D ASEP model with periodic boundary conditions describes a one-dimensional 
lattice of N sites, where each site i (1 < i < N) is either occupied by a particle or empty. 
During time interval dt, each particle has a probability dt of jumping to the adjacent site 
to its right if that site is empty. Our 2-D has a simplified dynamics since the time is now 
discretized (t„ = n). Each layer is a one dimensional lattice of sites, with periodic 
boundary conditions. The layers are labeled Lj. {k > 0). Each site is either occupied or 
empty and the density along each layer is set to be constant (= p). A shift of one half of 
the grid spacing is applied between consecutive layers, so that for k even i is an integer 
and for k odd, i is an half integer (see figure p. The 2-D lattice is composed of an infinite 
number of layers occupying the half space y > 0. At time t„, each particle has a probability 
Pi,k{tn) of hoping to the right. The quantity Pi^k{tn) is determined by the dynamics of the 
four nearest neighbors of the site {i, k): two in the row at time tn and two in the row 
Lk+i at time t„_i. Pi^kitn) is simply proportional to the number of these neighbors that 
are moving with a proportionality coefficient a. In addition, the exclusion principle imposes 
that the particle cannot jump to an occupied site. However, if the particle on site {i + 1, k) 
is jumping to its right at time t„ then the particle on site (z, k) is allowed to move. If Qi^k{t) 
is the characteristic function of the motion for the site {i, k) at time t {Qi,k{t) = 1 if there is 
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a particle at time t that is jumping from site {i, k), and is zero otherwise) the equation for 
Pi,k{tn) reads: 

Pi,k{tn) = Ot{Qi~l/2,k-l{tn) + Q i+l/2,k-l{tn) + Qi-l/2,fc+l (^n-l) + Qj+1/2,A:+1 (^n-l) ) (1) 

For consistency, if the right hand side of the formula (|I|) is bigger than 1 then we define 
Pi,k{tn) = 1- Eventually, our boundary condition will represent a moving wall situated at 
k = —1: all the sites of L_i are filled and are moving at each time step. No boundary 
condition is required for A; = oo. The velocity Vk{tn) is defined as the probability that a 
particle located on row L^. moves at time Then, the mean value of the velocity on row 
Lfc is denoted V{k) and defines the velocity profile. 

The only control parameter of the dynamics is a. Experimental observations suggest 
that p should be taken close to 1. It appears from the simulations that the dependence 
on p is trivial, and we consider results for p = 0.9. On the other hand, a reveals how 
much a moving particle pushes on its neighbors. Therefore, it has to be a (non trivial a 
priori) function of the material properties such as the friction, the shape of the grains, their 
roughness, etc... 

The model does not allow exchange of particles from one row to another (along the 
y direction). This strong constraint is in opposition with the experimental observations, 
where the density profile seems to reach a steady state where the exchanges between rows 
just balance. We assume in fact in this model that this interchange of particles is not 
relevant compared to the friction effects. We also checked that by imposing a reasonable 
density profile from the moving boundary to the bulk, the qualitative results of the model 
are not affected. 

III. NUMERICAL RESULTS. 

Figure shows different velocity profiles for p = 0.9 and a increasing from 0.2 to 0.65. 
For a > 1/2, an abrupt change in the velocity profile is observed. In fact, we can expand 
the logarithm of the velocity in a Taylor serie: 
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ln{V{k)) = a + b-k + c-k'^ + d-k^... (2) 

We remark that the three first terms on the right hand side of (0) give a qualitatively 
good approximation of the velocity profiles. It correspond to a Gaussian fit of the profile. 
However, the ratio b/c indicates the typical width for which the quadratic term becomes of 
the same order as the linear one. For a < 0.5 this ratio is of order of hundreds, i.e. much 
larger than the shear width and the dynamics can be considered almost purely exponential. 
On the other hand, the ratio b/c becomes of the order of 1 when a crosses the critical value 
ac = 0.5, so that one can approximate the velocity profiles for a > ac with a Gaussian 
centered near the row k = 0. Such property appears clearly on the insert of figure @, 
where the velocity profile for a = 0.65 as a function of k'^ is shown. Notice that for small 
k [k"^ < 500), the logarithm of the velocity is almost linear in k'^, so that one can consider 
that the velocity profile is Gaussian at least near the wall. 

The instantaneous velocity Vk{t) shows also different behaviors wether a is larger or 
smaller than ac. Thus, for a > 0.5, transitory dynamics occur for small time (t„ < 100) 
until the velocity reaches a stationary behavior. Such transitory states cannot be seen for 
a smaller than 0.5. 

In order to investigate more carefully the transition occurring at a = ac, as well as 
the velocity profiles, a simplified version of the model is introduced. It exhibits the same 
properties as the model explained above, and can be more easily studied analytically. It 
consists of neglecting the effect of the row L^+i on the dynamics in row L^. The model 
breaks the symmetry along the y direction. However, experimentally, the sheared flows 
exhibit a spontaneous symmetry breaking as well. A transition from exponential to Gaussian 
velocity profile occurs in the same way for this simplified model at a = ac. Notice that ac 
corresponds to the value of a at which the probability becomes 1 to move if the two neighbors 
above are moving (the exclusion condition still holds). 

We define P{n) as the probability that, given a hole, the next hole on its right is located 
after n filled sites. P{n) gives the probability distribution function (PDF) for the size of the 
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clusters (if two consecutive sites are empty, it is considered as a cluster of size 0). A priori, 
the function P{n) depends on p, a and the row number k. However, if we consider that the 
particles are placed randomly on the sites with a density p (as do the initial conditions), 
one obtains the Poisson distribution Po{n) = (1 — p)p"- As shown on figure (^, for a < ac, 
P{n) corresponds almost exactly to Po{n) for any k and a. Therefore, at each time step, the 
particles are distributed on the sites as if they were randomly placed with density p. On the 
other hand, for a > ac, the function P{n) differs from Po{n) in that the small size clusters 
are more frequent, while large clusters follow a Poisson-like law. In this case, P{n) does not 
vary as a increases for a given row, but changes as the row number increases: the bigger k, 
the closer P{n) approaches Po(^)- 

When a approaches ac the PDF differs slightly from Po{n), so for A; = 1 we can expand 
P{n) in the form (see insert of figure |^): 

P{n) = a6n,0 + (1 - a)Pe//(l - Peff), 

with pcff = p/(l — a(l — p)) (mean density has to be p). Figure (^) shows the evolution of a 
as a function of a for p = 0.9. For a > ac, a is a constant since a = Oc = 0.25 ± 0.01 which 
corresponds to Pe// ~ 0.92. However, for a < etc we observe that for a ~ ttc- 

ac — a (x {ac — aY 

and we found z/ ~ 2/3. Therefore, we identify the behavior of the system when a reaches 
Ol-Q clS 3i phase transition. Notice that when a increases above ac, the PDF shows a more 
complex behavior since not only the value of -P(O) is disturbed, but also a larger range of 
small clusters size n (see figure ||). 

IV. MEAN FIELD THEORY 

As the correlations can be neglected for a < ac, a mean field approach is appliable. 
Defining P^^\n) as the probability that n successive particles are moving on the row 
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(knowing that an empty site is before such a cluster), P^^\n) can be computed exactly: 
P(°)(ra) = (2ap)"(l - 2ap) and by noticing that 1/(0) = ^^J2nP^^'^{n), we obtain: 

-(°)^^ ' 

It is remarkable that the numerical simulations and this mean-field solution agree within 
an error of less than one percent. Also, one can write the constitutive relation between P*^°^ 
and P(^): 

oo 

pW{n)= W{1 ^n)P^'>\l) (3) 

l=n-l 

where W{1 ^ n) is the probability that if a cluster of size / at row /c = is moving, then a 
cluster of size n at row A; = 1 is moving. It follows that: 

W{l^n) = {l-p){2ap)"{{l-n + l)-{l-n)2ap) for l>n 

W{n ^ n) = (1 - p)ap{2ap)''-\2 - ap); W{n - 1 ^ n) = (1 - p){apy{2ap)''~^ 
and we obtain for V{1): 



P ^1 ' ' (l-2ap)2(l + 2ap) 

Again, this mean-field solution is in good agreement with the numerical results, although 

not as accurate than for V{0). Unfortunately, the next order step, for obtaining the analytical 
solutions of P^^^ and then V{2) are much more complicate. However, one can consider that 
for the exponential behavior, the knowledge of V{0) and V{1) is sufficient. Then, we can 
compare the numerical exponential profile with the exponential law predicted by the mean- 
field. We define r{a,p) such that for a < ac, the velocity profile obtained numerically is fit 
by: V{k) = V{0){r{a,p))'' 

Figure (^) shows log{r) as function of a for p = 0.9, compared with the mean-field r.^/ 
solution taken as: 



r 



mf 



V{1) _ ap{l - p)(2 + ap{l - 2ap) 
7(0) ~ 1 - (2ap)2 



Thus, the mean-field approximation, deducing the exponential law from the ratio be- 
tween the two first velocities, gives a quantitative good approximation of the exponential 
decay for a < ac- 

But, for a > the mean-field approach fails and we are not able to show analytical 
results. Although we were able to quantify the coorelations for a ~ a^, numerical simulations 
show correlations along the y direction as well. However, we believe that the phase transition 
exhibited by the model might have interesting features in granular fiows. Particularly, it 
might be interesting to perform an experiment where the friction coefficient of the granular 
materials would change. Also, the model shows a phase transition between exponential 
and Gaussain behavior as a increases, while experimentally, the different behaviors occur 
between 2 and 3 spatial dimensions. In fact, one can imagine, for example, that the 3-D 
experiment can be linked to this 2-D ASEP model with an effective friction coefficient Osd- 
Then, if one consider that in 3-D each particle has more neighbors that can push it, it is 
plausible to assume that a^o > and possibly a^n > aca2D under certain conditions. 
Also, we would like to point out other applications of this model such as non-newtonian 
fiows or molecular frictions. 

Finally, notice that he model is based mainly on probabilistic properties of granular fiows. 
Other stochastic approaches have already been proposed in this context although 



in granular materials there is no justification such as thermal fiuctuations for stochastic 
processes. However, the shape of the grains, and therefore the contact network in granular 
materials can be considered as random variables (this has been argued for the chain forces 
in static bead pile Additionally, the motion of the particles in shear fiows changes the 

configuration of the contact network of the system. 
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FIG. 1. 2-D lattice with sites (circle) at time the white circles represent the empty sites or 
holes, while the black circles represent particles. The motion of the particle located at the position 
{i, k) is determined as indicated by the motion of its four closest neighbors on rows L^+i and 
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FIG. 2. Velocity profile for a density p = 0.9 for different values of the control parameter a. 
a increases from 0.2 to 0.65 in increments of 0.05 as the curves are plotted from the left to the 
right. For a smaller than one half the profile is almost exponential and it is far from exponential 
for a > 0.5. The insert shows the velocity profile for a = 0.65 as a function of k'^. It shows that the 
profile is better approximated by a Gaussian law, particularly for the low values of k. The number 
of sites per row is 20000 and the average is computed over 4000 time iterations for 5 different initial 
conditions for each a. 
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FIG. 3. The probability distribution function of the cluster size for different values of a and 
k, and for p = 0.9. The straight line corresponds to Po{n) and is P{n) for a = 0.45 and k = 1. 
The two others PDF shown are for value of a larger than ac- The one that is maximum for n = 
is for a = 0.5 and k = 1 while the nest highest is for a = 0.6 and k = 4. Notice that they show 
Poisson-likc tails at large n, with different slopes. Insert: for a = 0.49 and k = 1, the PDF differs 
from Po{n) only at n = 0. 
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FIG. 4. The density of holes a as a function of a near the phase transition {a = Uc = 0.5). a 
is approximately constant for a > ac and exhibits a critical exponent of 2/3 for a < ac- 
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FIG. 5. log{r{p,a)) obtained by numerical simulation (circles) for p = 0.9 compared to the 
mean-field approximation logiV (1) /V {Q)) (line). log{r) is computed by a Gaussian fit of the 
numerical results. The error bars have been evaluated by comparing it with an exponential fit. 
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